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Abstract 

Coherent control of harmonic generation was studied theoretically. An optimal control theory 
was employed to find the driving field where restrictions were put on the frequency band. Addi- 
tional restrictions were added to suppress undesired outcomes such as ionization and dissociation. 
The method was formulated in the frequency domain where a relaxation updated procedure was 
' (— i ' employed. The method was tested on several examples demonstrating the ability to generate higher 

frequencies than the restricted frequency band of a driving field. 
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I. INTRODUCTION 



When high intensity light is irradiated on an atomic or molecular gas, higher frequencies 
are generated [1, 2]. This phenomenon is known as high harmonic generation [3-5]. The 
current approach to attosecond pulse generation [6] is based on this phenomenon. Significant 
effort has therefore been devoted to optimizing the process of harmonic generation [5, 7, 8]. 
Optimizing by control of the phase and amplitude of the incident light suggests coherent 
control [9]. 

The present paper address theoretically the issue of an optimal control strategy for har- 
monic generation. The target of optimization is the system's dipole operator. To be a 
source of radiation the acceleration of the dipole operator should oscillate in the frequency 
which is much higher than the driving field frequency. The idea is to exploit the significant 
theoretical development in optimal control theory OCT [10-14]. The hope is that the opti- 
mized pulses will unravel specific mechanisms of harmonic generation. OCT for harmonic 
generation poses two significant challenges: 

• A constraint on the bandwidth of the incident control pulse has to be imposed; 

• The target is designated in the frequency domain while OCT is typically formulated 
in the time domain. 

Several suggestions for dealing with the first issue appear in the literature, in a more general 
context [13-20]. It has been attempted to deal with the second issue by the means of the 
general OCT formulation of time-dependent targets [13, 21]. However, no results from this 
approach have been reported. 

In the present study, the two challenges are overcome by the formulation the control 
problem in the frequency domain, replacing the formulation in the time domain. In this 
approach, the frequency requirements of the problem are expressed in a natural and direct 
way. A simple and effective optimization procedure, suitable for the new formulation, is 
suggested. 

II. OCT OF TIME DEPENDENT PROBLEMS 

We first review the formulation of optimal control theory for time dependent problems. 
Let us denote the time-dependent state of the system by \ip(t)), the drift Hamiltonian by H , 
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and the driving field by e(t). The dynamics of the system is governed by the time- dependent 
Schrodinger equation, under a given initial condition: 

am)) 



dt 



\m) = i^ > 



(i) 



where H(t) = H — p>e{t). (Atomic units are used throughout, so we set: h = 1.) The 
optimization functional for time-dependent targets becomes (see [13, 21-23]): 



J — J max ~\~ Jbound ~\~ Jpenal ~\~ J a 
Jmax = J W(t)(ll>{t)\6(t)\ll>{t)) dt 

J bound = k(^(T)\6{T)^{T) 

Jpenal = ~Oi I (t) dt 

JO 

J ran 



w(t) > 0, 



w(t) dt 



K>0 

a > 



T 



e/ Q {X(t) 







Dt 



+ iH(t) 



7p(t) ) dt 



(2) 
(3) 
(4) 
(5) 

(6) 



where 0(t) is the time-dependent target operator, and T is the final time. J ma x represents 
the target to be maximized. Jbound is a boundary term. The inclusion of this term prevents 
boundary problems (see [23, Section 2.2]). J pena i is a penalty term on the intensity of 
e(t)- J con represents the constraint on the dynamics of the system — the time-dependent 
Schrodinger equation. 

The resulting Euler-Lagrange equations become: 

d\m) 



dt 
dixit)) 
dt 



-m(t) \m) , 

-m(t) \ x (t))-w(t)6(t) \m) 



\m) = i^o) 

\X(T)} = k6(T) \rf,{T)) 



(7) 



ft(t) = Ho-/2e(t) 



a 



(9) 



This formulation is suitable for targets that are well defined in the time domain. For prob- 
lems with frequency requirements, this approach has to be modified. 

III. OCT OF HARMONIC GENERATION 



The harmonic generation problem may be divided into two distinct parts: 
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1. The driving field spectrum has to be restricted to the frequency range available from 
the source; 

2. The intensity of the emission of the system in the desired portion of the spectrum has 
to be maximized. 

These two parts will be treated separately in the next two subsections. 
A. Restriction of the driving field spectrum 

The task of restricting the driving field spectrum has considerable importance in OCT; 
the reason is that most of the computed fields turn out to be too oscillatory to be produced 
experimentally. This problem may be overcome by limiting the spectrum of the field to 
sufficiently low frequencies. 

Several approaches for achieving this goal have been proposed. In [13], a spectral fil- 
tration of the driving field is performed in every iteration of the Krotov algorithm. This 
approach leads to a non-monotonic convergence of the optimization procedure. In [14], a 
two-dimensional penalty term is introduced, in order to control the spectral properties of 
the driving field. This approach might lead to numerical instabilities, or non-monotonic 
convergence of the optimization procedure (see [23, Section 3.2.1]). In [15], the problem 
of optimization of a general driving field function is replaced by the optimization of the 
coefficients of a list of frequency terms. 

In the present approach, the restriction on the field spectrum is achieved by placing a 
penalty function on the undesirable frequency components of the field. The regular J pena i 
from Eq. (5) is replaced by a penalty term formulated in the frequency domain. In [23, 
Section 3.1.2], it is shown that there is a close relationship between this formulation to the 
methods presented in [13-15]. 

The cosine transform is employed as a spectral tool. Other spectral transforms (i. e. the 
Fourier transform or the sine transform) could be used as well. For a typical signal, a cosine 
series is known to converge faster than a Fourier series or a sine series (see [23, Section 3.1.1]). 

The operation of the cosine transform on an arbitrary function g(t) is denoted by the 
symbol C, and the transformed function by g(oS): 




(10) 
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The inverse cosine transform will be denoted by C 



-l. 



C- l [g(uj)} = \l- I g(u) cos(cot) du = g(t) (11) 



The driving field in the frequency domain is defined by a finite time cosine transform: 

e(u) = \[^- [ e(t) cos(wt) dt (12) 



The maximal cutoff driving frequency is denoted by Q. The penalty term from Eq. (5) 
is modified to: 

Jpenal = -(y -—-e 2 (u)du «>0 (13) 

JO fe{u) 

where f e (u>) is an adjustable function, which satisfies the conditions: 

f e (u)du = l, / e (w)>0 (14) 

a determines the cost of large fields, Cf. Eq. (5). f e {w) is chosen to have small values 
for undesirable frequencies, and regular values for the allowed frequency region. It may be 
interpreted as a filter function, as can be seen in Section III C. An additional envelope shape 
can be forced on the profile of the driving field spectrum by choosing an appropriate filter 
function. A complete filtration of undesirable frequencies is achieved in the limit f t (uj) — > 
for undesirable oo values. For practical purposes, f e (oj) may be set to for these values. 



B. Optimization functional for harmonic generation 

The field emitted by the system consists of the frequencies contained in the spectrum of 
the dipole expectation. In order to maximize the emission in a desired frequency region, the 
amplitude of the dipole moment oscillations in this frequency region is maximized. Thus, 
the physical quantity of interest is the dipole moment expectation value: 

(/*)(*) = Mt)|AlV(*)> (15) 

The dipole spectrum becomes: 

Jfi(u>)=C[(fl)(t)] = yf £ (fi)(t)cos(u;t)dt (16) 
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To maximize the emission in the desired region of the spectrum, the following functional is 
chosen: 

1 f a 2 

Jrnax = ->^ f ^(u) (ft) (u) d(J A>0 (17) 

z Jo 

where satisfies the conditions: 

f,(u)du = l, /»>0 (18) 



o 

A is an adjustable coefficient, which determines the relative importance of J m ax- A is redun- 
dant with a but from numerical stability considerations it is useful to vary it independently. 
fn(uj) is a filter function, which has pronounced values in the frequency region of interest. 
Eq. (17) can be generalized to an arbitrary Hermitian operator O: 

n 



o 



J max = -Xl /oH(O) (u)du; A>0 (19) 



0)M =C 

n 



O) (t)\ (20) 
fo(u)du = l, /oH>0 (21) 



o 



One possible application of this generalization is in the case when it is desired to maximize 
emission with polarization other than that of the control field. For instance, if the control 
field is x polarized, and we require emission of y polarized field, p, in H(t) is set to be £i x , 
and O = p, y . 

The full maximization functional for the harmonic generation problem becomes: 

J = Jmax Jpenal Jcon (22) 

1 f n Try 2 
J max = -\J f (u)(p) (cu)du (23) 

Jpenal = -a —- 6 2 (uj) du (24) 







J con = -2Re ( X (t) 



T 



i/j(t) ) dt (25) 



C. The Euler-Lagrange equations for harmonic generation 

We choose the functional derivative of the objective in the frequency domain: 
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The resulting Euler-Lagrange equations become: 

d\m) 



dt 
dixit)) 
dt 



-m(t) m)) , 

-iU(t) \ X (t)) - XC' 1 



fo(u)(6)(u 



o\m) 



\m) = 1^0) 

\X(T)} = 



(27) 
(28) 



where H(£) = H — p,e(t) and 

e(u) = Mu)C[r!(t)] t 



rj(t) 



lm( X {t) \fi\m) 



a 



e(t)=C- 1 [e( U )]=C- 1 {fMC[v{t)]} 



(29) 
(30) 



Note that the expression for rj(t) is the same as that for e(t) in Eq. (9). e(t) in Eq. (30) 
can be interpreted as the filtered field from the regular control problems, where f € {oj) plays 
the role of a filter function. A comparison between Eq. (8) and Eq. (28) leads to a similar 
interpretation of the inhomogeneous term in Eq. (28) (Cf. [23, Section 3.3.1]). 

It is convenient to avoid normalizing f e (oj) and fo{w), and substitute them with: 



a 



fo(u) = Xf (u) 



(31) 
(32) 



D. Optional modifications for the optimization problem 



1. Prevention of dissociation 

Typically when strong driving fields are employed the system dissociates or ionizes. This 
phenomena can be avoided by restricting the system to localize in an "allowed" subspace 
of the Hilbert space. For example, eliminating the access to all eigenstates with energies 
above a threshold energy. Another option is to restrict the state vector to regions of space 
far from the threshold of the potential well. A similar idea n suppressing dissociation is to 
restrict the dynamics to the allowed momentum values. 

In order to restrict the system to the allowed subspace, two modifications in the maxi- 
mization functional J are employed: 

1- J max is modified to include contribution only from the allowed states; 

2. A penalty term on the forbidden states is added to J. 
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The first modification is achieved by the replacement of the expectation ( O ) (t) in 
Eq. (20) by the expression: 

'P a if>(t) O P a ^(t) 



where P a is the projection operator onto the allowed subspace. For instance, if all energies 
above the threshold El are restricted, then: 

L 

Pa = ^ ^ ( 33 ) 
n=0 

If the system is restricted in the x space to remain in the interval [x m i n , x max ], then: 

r^max 

P a = \x) (x\ dx (34) 



If a smooth filtration of states is desired, P a may be generalized to a weighted projection 
operator, that will be denoted as P*. For instance, P a from Eq. (33) is modified to: 

JV-l 

P s a = J2s n MM 0<s n <l (35) 

n=0 

where s n decreases gradually from 1 to near the threshold. P a from Eq. (34) is modified 
to: 

K = J s O) \ x ) ( x \ dx = s(lCj < s(x) < 1 (36) 

where s(x) decays to near the boundaries of the allowed x interval. 
The resulting modified J max becomes: 

i rtt —, —2 





6 = p s op s 



Jmax = -I /o(w)(O a ) (u)dlU (37) 



The second modification is the addition of the following penalty term to J (as suggested 
in [22]): 

J f orb = -lJ (^)|P/|V>(0) dt 7 >0 (38) 

where P/ is the projection onto the forbidden subspace, and 7 is the penalty factor of the 
forbidden subspace. 
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It is possible to achieve a smooth filtration of states by using a state-dependent penalty 
factor. For instance, in the energy space Eq. (38) is generalized to: 



J 



forb 



p7 



^2 In \<f n ) (<Pn 



p7 



ij>(t)) dt 



N-l 



In > 



(39) 
(40) 



n=0 



Pj is a skewed projection. 7„ is the penalty factor of the state \<p n ). It should be for the 
allowed energy domain, and increase gradually with n near the threshold energy. In the x 
space, the skewed projection is: 



P] = / 7(2) \x) (x\ dx = 7 (X 



l{x) > 



(41) 



where 7(2) increases gradually in the forbidden x regions. A smooth penalization is recom- 
mended, in order to decrease the difficulty in the optimization process. 

When the Hilbert state is very large, it becomes impractical to compute all the eigenstates. 
Nevertheless, a restriction of the allowed subspace in the energy space is still possible. s n 
and 7 n should be defined as functions of the energy, i. e. : 



s(E n 



In = l{E n ) 



Then we have: 

K = 8 (Ho) PJ = 7 (Ho) (42) 

s (jloj \i>{t)) an d 7 (-^-o) IV ; (^)) can be approximated using standard methods. 
When P^ = i and P] = 6, J reduces to Eq. (22). 

After inserting these changes in J, the equation for \x(t)) Eq. (28) is modified in the 
following way: 

d\x(t)) 



dt 



-m) \x(t)) - c- 1 



fo(u)(O a )(u 



o a 



pt s> \m) 



(43) 



2. Prevention of boundary effects 



In some applications undesirable boundary effects at the final time interval should be 
suppressed (see [23, Section 3.3.3]). This is achieved by adding the penalty term: 

'd(6)(Tf 2 



Ji 



bound 



1 

■-K 

2 



dt 



k > 



(44) 
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k is an adjustable parameter. When k = 0, J reduces to Eq. (22). 



In the special case that p,, O 



6, the insertion of Jbound results in a relatively simple 
modification of the Euler-Lagrange equations. The natural boundary condition for in 
Eq. (28) is replaced by: 



\X(T)) = k 



H ,O 



(T) 



Hn,0 



(45) 



The general case is more complex, and will not be discussed here. 

The derivation of the Euler-Lagrange equations is presented in [23, Appendix A]. 



E. Optimization procedure 

For current optimization posed both in time and frequency the well established opti- 
mization procedures do not converge monotonically or converge very slowly. We therefore 
employed a more direct relaxation method to update the driving field from iteration to itera- 
tion. In the present context the method consists of the following update rule for the driving 
field: 

e new {u) =Ke EL {u) + {l-K)€° ld {u) < K < 1 (46) 



= f £ (uj)C 



-lm( X (t) 



(47) 



e(u})=e old (ui) 

The updated field is a mixture of the previous field, and the field computed from the Euler- 
Lagrange equation Eq. (29), using the previous field for the computation of \x(t)) an d \ip(t)). 
K is the mixing parameter, which determines the weights of the two fields. The value of K 
is decreased when the optimization process progresses. 

The scheme for the implementation of the relaxation method becomes: 

1. Guess a driving field spectrum e^(co). 

2. Set: e(°)(t) = C- 1 [e(°)(u;)]. 

3. Guess an initial value for K. 

4. Propagate \^°\t)) forward from t = to t = T according to Eq. (27), with e^°\t). 

5. Calculate with |^ (0) (t)) and e^(uj). 
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6. (k = 0) 

7. Repeat the following steps until convergence: 

(a) Set \x (k) (T)) according to Eq. (45), using \ip (k) (T)) and e (fe) (T). 

(b) Propagate backward from t = T to t = according to Eq. (43), with 
e( fe )(t). 

(c) Do the following steps, and repeat while J tnal < j( fc ); 

i. Set a new field, using Eq. (46): 



e*™>) = Kf e (u)C 



lm(x ik) (t) 



^ k \t) 



[1-K)e {k \u 



e trml (t) =C- 1 [e tria \u)] 



ii. Propagate m trml (t)) forward from t = to t — T according to Eq. (27), with 
e tria \t). 

iii. Calculate J trial with \ip trial {t)) and e trial (u). 

iv. If J trial < J( fc ), then set: K = K/2 

(d) Update all the variables: 

^ k+1 \u) = e trial (u) e {k+1) (t) = e trial (t) |^ (fc+1) (^)> = \^ tria \t)) J (k+1) = J trial 

(e) (k = k + 1) 

In [23, Section 3.2.3] it is shown that the relaxation method, in the context of quantum- 
OCT problems, is equivalent to a second order gradient method (quasi-Newton method). 

IV. APPLICATION 

The new method is demonstrated in four simple harmonic generation examples. The 
propagator for the Schrodinger equation is based on a new, efficient and highly accurate 
algorithm, [24]. 

The convergence condition of the optimization procedure is: 



I ^ new ^ old I 



< t (48) 
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0.5 1 1.5 2 

u> (a.u.) 



FIG. 1: J^uj) of the TLS problem 

where e is the discrete vector of frequency values on an equidistant u grid and r is a tolerance 
parameter. More numerical details may be found in [23, Appendix B]. 

The important details of the problems and the computational process are presented in 
the tables. Atomic units are used throughout. The notations in the tables are described 
in Table I. O(x) denotes the Heaviside step function. The initial state in all problems is 
chosen to be the ground state, denoted as |<^o)- 

A. Two level system 

The first problem is tripling the driving frequency by a two level system (TLS). The 
unperturbed Hamiltonian is: 

(49) 

The dipole moment operator is chosen to be the x Pauli matrix: 

(50) 

The driving field is restricted to be centered around u = l a . u ., by a "hat" filter function 
(Cf. Fig. 1). We require maximization of the emission in the region of the characteristic 
frequency of the system, = 3 a . u .. A Gaussian function is used for f^(u) (see Fig. 2). 
The details of the problem are summarised in Table II. 



Hn 



1 
4 



1 

1 
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Notation 


Description 


H 


unperturbed Hamiltonian 


A 


dipole moment operator 




initial state vector 


T 


final time 




scaled filter function of the driving field 


IM 


scaled filter function of the dipole moment expectation value 




initial guess of the field 


L 


index of the maximal allowed eigenstate 


7n 


penalty factor of the state \ip n ) 


s(x) 


projection function onto allowed x regions 


7(x) 


penalty function on forbidden x regions 


*Q 


initial guess of K, for the relaxation method 


x domain 


domain of the x grid 




number of equidistant points in the x grid 


T 


tolerance parameter of the optimization process 



TABLE I: Description of the notations in the tables 

The optimization process converges rapidly to a solution. The convergence curve is shown 
in Fig. 3. 

The resulting spectra of the driving field and the dipole moment expectation value are 
shown in Fig. 4. e(u) is shown to be successfully restricted to the desired portion of the 
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to (a.u.) 

FIG. 2: f^u) of the TLS problem 



H 


ri 01 
L 4 J 


A 




IV>o> 


[h] 


T 


100 


AH 


20sech[20(u; - l) 4 ] 




exp[-10(£j - 3) 2 ] 




sech[20(a; - l) 4 ] 




0.5 


r 


10" 3 



TABLE II: The details of the TLS problem 

spectrum. The "hat" envelope shape is apparent. (p)(u)) mainly consists of a large peak at 
cu^o, as required. 

B. Eleven level system 

The second problem is of an eleven level system (11LS). The problem is designed for 
harmonic generation by a resonance mediated absorption mechanism (for example, see [25]). 
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1 2 3 4 5 

number of iterations 



FIG. 3: The convergence curve of the TLS problem 



8 
6 

_ 4 
3 

d 

X 2 



-2 



12 3 4 

uj (a.u.) 



FIG. 4: The spectra of the driving field (red) and the dipole moment expectation spectra (blue), 
for the TLS problem 
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The unperturbed Hamiltonian is: 



H r 



2.1 







3.9 

5 

6.1 

7 

8.1 



9.9 
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The dipole moment operator is: 



1 

1 1 

1 1 
1 1 



10 1 
1 1 
^1 



1 



1 1 
1 1 
1 1 
1 



(51) 



(52) 



This jl couples between neighbouring eigenstates, and between the outer eigenstates, \ifo) 
and Iv^io)- The driving field is restricted not to exceed the region of the resonance frequencies 
of the neighbouring levels, u) n +i,n = la.«.- We require an emission in the neighbourhood of 
the Bohr frequency of the outer energy levels, 6^10,0 = 10 a .u.- fe(u) and / M (w) are chosen to 
be rectangular functions. 

The details of the problem are summarized in Table III. 
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H 


Eq. (51) 




Eq. (52) 


1 / \ 


1 \ 

m) 


T 


100 




50 6(1.3 — ui) 




Q(io - 9.9) 9(10.1 - oj) 




fc)(l.o — w) 


Ki 


1 


T 


10~ 3 



TABLE III: The details of the 11LS problem 



8 

7 
6 
5 

~i 4 

3 
2 
1 


50 100 150 

number of iterations 

FIG. 5: The convergence curve of the 11LS problem 

The convergence curve is shown in Fig. 5. The resulting e(oS) and {p)(u) are shown in 
Fig. 6. The new method is shown again to be quite efficient in maximizing the emission in 
the required region. 
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10 
5 

° i 

-5,1 

3 

-10 
-15 

2 4 6 8 10 12 

lo (a.u.) 

FIG. 6: The spectra of the driving field (red) and the dipole moment expectation value (blue), for 
the 11LS problem 

C. Anharmonic oscillator — the HC1 molecule 

In this example an anharmonic oscillator is used to double the incoming frequency. The 
oscillator is chosen is an approximation of the H 35 C1 molecule (see [23, Section 4.3.1] for 
more details). As in the previous problem, the intended mechanism is of harmonic generation 
by resonance mediated absorption. 

The coordinate of the one- dimensional oscillator is the displacement of the inter-nuclei 
distance, th-ci, from the bottom of the well (r*): 

x = r H -ci - r* (53) 

The approximated potential function, V(x), and dipole moment function, fi(x), are shown 
in Fig. 7. 

e(oj) is restricted not to exceed much the characteristic frequency of the bottom of the 
well: 

u = 1.35 • 10" 2 a . u . 

We require maximization of the emission in the neighbourhood of the second harmonic: 

w 2 ,o = E 2 -E = 2.54 ■ 10~\. u . 
f € {oj) and fn(u)) are chosen to be rectangular functions. 
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1 2 3 

x (a.u.) 

FIG. 7: The approximated potential (blue) and dipole function (red) curves, for the HC1 molecule 




20 40 60 80 100 120 140 

number of iterations 



FIG. 8: The convergence curve of the HC1 problem 

In order to prevent dissociation of the molecule, the energies above the dissociation thresh- 
old were restricted (see Section IIID1). 

The details of the problem are summarised in Table IV. 

The convergence curve is shown in Fig. 8. The resulting e(u) and (fi)(oj) are shown in 
Fig. 9. (/t)(w) mainly consists of a large linear response to the driving field, as could be 
expected. However, there is also a significant non-linear response in the region of the second 
harmonic frequency, as required. 

More information on the results of the first three examples, and a detailed discussion, 
may be found in [23, Chapter 4]. 
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H 


2 f 7 2 g5 + 0-171 [exp ( 0.975 x) i] 2 




(0.19309 x) x 

jl-Re tanh ^(0.17069 + 0.056854 i) (± - 0.10630 f) 1 ' 8977 ^ 1 


l</>0> 




T 


10 4 




2500 9(0.015 - oj) 




100 e(w - 0.025) 6(0.027 - w) 


L 


19 


In 


f n<19 
\ (n-19) 2 n>19 




8(0.015 -w) 




1 


x domain 


[-0.69407, 3.51178) 




32 


r 


10^ 3 



TABLE IV: The details of the HC1 problem 
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300 




200 



100 2 



4r 



o 



3 



-100 



-200 



-300. 







0.005 



0.01 



1 0.015 

lu (a.u.) 



0.02 



0.025 



FIG. 9: The spectra of the driving field (red) and the dipole expectation spectrum (blue), for the 
HC1 problem. Two significant peaks appear: a large linear response of the dipole to the driving 
field, and a significant non-linear response in the neighbourhood of the second harmonic. Notice 
the low frequency component of the driving field which changes the static part of the Hamiltonian. 

D. One dimensional particle in a truncated Coulomb potential 

The last example demonstrates the application of the new method in a system with 
stronger non-linearity. A driven electron in a Coulomb potential is the system studied. The 
model has been extensively studied in the context of harmonic generation [4, 5]. Typically 
for strong driving fields a comb of odd frequencies is generated up to a cutoff. In the present 
example the target is the emission of a single high harmonic of the driving frequency This 
target has similarities to the experiment in JILA [26] where the emission of the 27 harmonic 
was enhanced relative to its neighbours using a genetic algorithm optimization. 

Our model consists of a particle of unit mass and charge placed in a truncated Coulomb 
potential constrained to one dimension (see Fig. 10): 



The driving field is restricted to frequencies which are much lower than the resonance 
frequencies of the system. e(u) is restricted not to exceed the region of u — 0.07 a . u .. The 
shape of f e (u) (see Fig. 11) induces a smooth filtration of higher frequencies. We require 
maximization of the emission in the region of the one of the Bohr frequencies of the system, 




(54) 



The dipole operator is /t = X. 
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FIG. 10: The truncated Coulomb potential (Eq. (54), black), and the potential energy of the 
system under the influence of a strong field (dashed red). 

x 10 4 




0.02 0.04 0.06 0.08 0.1 0.12 

uj (a.u.) 



FIG. 11: / e (w) of the truncated Coulomb potential problem 

^5,o = 0.624 a . u .. For the filter f^(u) a rectangular function is employed. 

The edges of the x grid are restricted using the method presented in Section HID 1. s(x) 
and 7(x) are shown in Fig. 12. 

The details of the problem are summarised in Table V. 

The convergence curve is shown in Fig. 13. 

The resulting field e(u) and dipole (£l)(uj) are shown in Fig. 14. The largest peak of (£i){oj) 
is located near the fundamental frequency of the system, = 0.395 a . u .. The response in 
the desired frequency is marked on the figure. The method is shown to be effective also in 
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FIG. 12: s(x) and 7(2;) of the truncated Coulomb potential problem 
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FIG. 13: The convergence curve of the truncated Coulomb potential problem 

the production of frequencies considerably higher than those of the driving field. 

The final solution obtained up to the chosen r is not a converged solution, as can be 
deduced from the convergence curve shape at the end of the optimization process (see 
Fig. 13). If the optimization process is carried on to smaller r, contributions to (£l)(oj) 
from interactions with the boundaries of the x grid begin to be significant. Nevertheless, 
a converged solution which do not involve such undesirable effects can be obtained. This 
is achieved by the increment of the j(x) values during the optimization process, before the 
spurious effects appear. 

Maximizing the response of an arbitrary frequency, which is not one of the Bohr fre- 
quencies, was also achieved. However, the resulting amplitude of the response was found 
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TABLE V: The details of the truncated Coulomb potential problem 

to be considerably smaller than if a Bohr frequency was selected. Larger response was ob- 
tained when partial ionization was allowed. Technically this requires absorbing boundary 
conditions. This topic is still under investigation, and therefore is not presented. 

V. CONCLUSION 

Optimizing harmonic generation is one of the most difficult tasks in the context of quan- 
tum control. A major obstacle is that the target objective cannot be formulated in the time 
domain. Additional restriction have to be added to suppress ionization or dissociation. A 
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FIG. 14: The spectra of the driving field (red) and the dipole expectation spectrum (blue), for 
the truncated Coulomb potential problem; the response in the desired frequency, u^o = 0.624 a u . , 
is marked by a black ellipse. The largest peak of (p,)(u) is in the fundamental frequency, 
win = 0.395 a . u .. 



new theoretical method of calculation optimal control of harmonic generation was studied. 
The task was addressed by the means of OCT, using a frequency domain formulation. The 
relaxation method was used as the iterative optimization procedure. 

For low and intermediate control fields fast convergence was obtained when the emitted 
high harmonic fit a fundamental Bohr frequency of the system. Stronger fields can modify 
the system thus allowing harmonic emission in other frequencies. The difficulty we found in 
locating such solutions was that they competed with dissociation or ionization. A suggestion 
in this direction is the low frequency contribution of the control field in the HC1 example 
which modified the Bohr frequencies (see [23, Chapter 4.3.1] for more details). 

The new method provides a useful tool for theoretical investigation of harmonic gener- 
ation mechanisms (as was demonstrated in [23, Chapter 4]). We expect that the current 
approach will contribute to the understanding of harmonic generation processes, and that 
new harmonic generation mechanisms will be revealed. 
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